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We perform molecular dynamics (MD) simulations with phonon spectral analysis aiming at under- 
standing the two dimensional (2D) thermal transport in suspended and supported graphene. Within 
the framework of equilibrium MD simulations, we perform spectral energy density (SED) analysis to 
obtain the lifetime of individual phonon modes. The per-mode contribution to thermal conductivity 
is then calculated to obtain the lattice thermal conductivity in the temperature range 300-650 K. 
In contrast to prior studies, our results suggest that the contribution from out-of-plane acoustic (or 
ZA) branch to thermal conductivity is around 25-30% in suspended single-layer graphene (SLG) at 
room temperature. The thermal conductivity is found to reduce when SLG is put on amorphous 
SiCb substrate. Such reduction is attributed to the strengthened scattering in all phonon modes in 
the presence of the substrate. Among them, ZA modes are mostly affected with their contribution 
to thermal conductivity reduced to around 15%. As a result, thermal transport is dominated by 
in-plane acoustic phonon modes in supported SLG. 



I. INTRODUCTION 

Ever since its experimental isolation in 2004 1 , graphene 
has been the focus of extensive studies because of a num- 
ber of fascinating properties it exhibits 2 4 . Aside from 
high electronic mobility 5 , 2D thermal transport and su- 
perior thermal conductivity have also been demonstrated 
in graphene experimentally 4,6 . It was suggested that such 
2D transport can lead to infinitely large intrinsic thermal 
conductivity - —. Also, it leads to novel three dimensional 
(3D) to 2D crossover phenomena in thermal conductivity 
when the material dimensions shrink from corresponding 
bulk 10,11 . In addition, when extrinsic perturbations, for 
instance, corrugation, grain boundaries, substrate, etc., 
are present, additional phonon scattering channels can 
be opened and contribute to the thermal resistance in a 
manner distinguished from those typically found in other 
material systems. 

Despite extensive investigations on the electronic 
transport in graphene, many issues of thermal transport 
still remain unclear. For instance, reported thermal con- 
ductivity values of suspended SLG spread widely between 
1,800 and 5,300 W/m-K at room temperature 4 ^^ possi- 
bly due to the varying degrees of sample quality, render- 
ing it controversial in identifying 2D thermal transport 
features in SLG. In order to address these issues, many 
theoretical approaches have been proposed 13 17 . Among 
them, Nika et aA 7 - suggested negligible contribution from 
ZA phonons due to their low group velocity. Lindsay 
et aA 5 - suggested, on the contrary, that the majority heat 
is carried by ZA phonons. It is evident that agreement 
on the relative importance of individual phonon modes 
to the thermal conductivity in suspended graphene has 
not been reached, largely because these models used dif- 
ferent assumptions that cannot be verified directly with 
experiments. This issue extends to the understanding of 



the reduced thermal conductivity observed in supported 
graphene 6 - while few quantitative theoretical models 18 
without adjustable parameters have been developed to 
predict its thermal conductivity so far. 

In this work, we employ MD simulations in combina- 
tion with SED analysis to extract spectral phonon life- 
time, mean free path (MFP) and thermal conductivity for 
both suspended and supported SLG in the temperature 
range 300-650 K. Contrary to previous suggestions, we 
find that the out-of-plane phonons (ZA and ZO) couple 
to the in-plane phonons in suspended SLG with medium 
strength, and the contribution of ZA phonon to the to- 
tal thermal conductivity is neither negligible nor domi- 
nant. We also find that the contributions from all phonon 
modes to thermal conductivity are reduced in the pres- 
ence of substrate while ZA modes are mostly affected. 



II. METHODOLOGY 

A. Molecular dynamics simulations setup 

In Fig. [I] we present the MD simulation domain setup 
for both suspended and supported SLG. The dimensions 
of the shown domains are 4.4 x 4.3 x 1.6 nm 3 for length, 
width, and height, respectively. Periodic boundary con- 
dition (PBC) is applied in the in-plane directions (paral- 
lel to graphene plane) while no restriction is applied to 
SLG vibration in the out-of-plane direction. We choose 
silicon dioxide (Si02) as the substrate, which is usually 
used in experiments. A 2 nm thick block of amorphous 
Si02 which comply with the PBC of the graphene domain 
is prepared by following the heating-quenching recipe 
used by Ong et al 19 . Similar thickness of substrate has 
been successfully used previously for both electronic^ - — 
and thermal simulations 18,19 . Thicker Si02 block is not 
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FIG. 1: (Color online) Instantaneous geometry of both sus- 
pended and supported SLG. The SLG plane is parallel to the 
x-y plane while the SLG out-of-plane direction aligns with z. 

found to induce obvious difference in the in-plane ther- 
mal transport of supported graphene. To ensure stability 
and to approximate semi-infinite substrate, the bottom 
layer of amorphous Si02 is fixed. The graphene sheet is 
released from 2 A above the substrate and allowed to con- 
form to the Si02 surface freely. The geometries of both 
suspended and supported graphene are pre-optimized un- 
der constant pressure and temperature (NPT) for 6 ns 
to ensure zeroed pressure, equilibrated temperature and 
stablized interfacial structure before being used for indi- 
vidual simulation runs. 

The LAMMPS package is used to perform all MD 
simulations 2 - 3 . Optimized Tersoff (OPT) potentials 24 are 
adopted to model C-C interactions in SLG. Interactions 
within Si02 are modeled using Tersoff potentials param- 
eterized by Munetoh et a/—. C-Si and C-0 interactions 
are assumed to be of van der Waals (vdW) type 2 ^ and 
modeled using Lennard- Jones (LJ) potentials 

y(r^) = 4e[(f- ) 12 -(^) 6 ] (1) 

with parameters ec-Si = 8.909 meV, crc-Si — 3.326 A, 
ec-o = 3.442 meV and a C -o = 3.001 A 26 . The cutoff 
is chosen to be 2.7a. Here nj is the separation between 
atom i and j. The timesteps for suspended and sup- 
ported SLG simulations are chosen to be 0.8 and 0.2 fs, 
respectively, to ensure stability and resolution of all possi- 
ble vibrational frequencies. Starting from pre-optimized 
geometries, after initial equilibration in constant volume 
and temperature ensemble (NVT), the systems are run 
in the constant volume and energy ensemble (NVE), 
from which the atomic velocities are extracted and post- 
processed. The phonon properties are found to converge 
for simulation lengths of 3.2 and 1.0 ns in NVE ensem- 
ble for suspended and supported SLG, respectively. In 
this work, only the thermal transport in x direction (zig- 
zag) is considered. Five independent simulations runs are 
performed and averaged for each case study to minimize 
effects from statistical fluctuations. Since the tempera- 
tures at which we perform the MD simulations are below 
the Debye temperature of SLG (~ 1000-2300 K— ), quan- 



tum corrections (QC) outlined in the Appendix have been 
performed, which were successfully practiced in similar 
ways in previous work o 13 i 28 . 

B. Spectral energy density approach 

In most cases, the intrinsic phonon transport in crys- 
tals is dominated by the three-phonon Normal or Umk- 
lapp process. Therefore, considering only up to third- 
order anharmonic force constants is usually a good ap- 
proximation. However, in CNT's, neglecting higher-order 
anharmonic interactions is found to significantly over- 
estimate phonon lifetime 29 . Due to the close connec- 
tion between CNT and SLG, we expect similar overes- 
timation in SLG if a fully anharmonic simulation is not 
performed. To predict phonon dispersion relations and 
lifetime of SLG, we employ the phonon SED analysis 
which, in combination with MD simulations, naturally 
includes the fully temperature-dependent anharmonicity 
of the atomic interactions. Such analysis and similar for- 
mulations have been successful in predicting individual 
phonon properties for various material systems 3 -^"— , and 
hence will be used here. 

Although phonon lifetime prediction using only atomic 
velocity data was reported in recent SED works on 
CNT 29 , PbTe 36 , suspended graphen o 37 i 38 , and supported 
CNT 39 , it was later pointed out by Larkin et al that 
eigen-displacements should also be included in the for- 
mulation of SED functions^. A detailed derivation of 
the connection between SED function and phonon spec- 
tral properties was given in Ref4£. Here we present a 
simplified version to illustrate the methodology. 

In a harmonic system, phonons do not scatter. There- 
fore they have infinite lifetimes and exhibit delta-function 
type peaks in vibrational spectrum. In real material sys- 
tems, phonons scatter and thus have finite lifetimes due 
to the anharmonic interactions, leading to shift in the 
phonon spectrum and broadening of the phonon peaks. 
According to Ladd et aJ2ft, the normal mode amplitude 
can be written under single mode relaxation time approx- 
imation as 

S k ,„(t) = S k)i ,, e-'«,-^,,)*. (2 ) 

Here So is the magnitude of the vibration, uj a is the 
anharmonic angular phonon frequency, Y is the phonon 
linewidth, k is the wavevector, v is the index of phonon 
branches. Then Fourier transform of the time derivative 
of normal mode amplitude is 

1 f°° 

F[s^(t)] = -= / &,„,(>(-««<„ - r k)V ) 

V Z7T J-oo 

xe-^^-^e-^^dt. (3) 

Its norm square is: 

9 Si „ Mud J 2 + T 2 )/2tt 

\ F [s K At)}\ 2 = 7^_^) 2 + r 2 • (4) 
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This function is in the Lorentzian form with its full width 
at half maximum (FWHM) equals 2T. In addition, the 
phonon lifetime r is directly related to the linewidth as 30 : 



r = 1/2T. 

Therefore, when SED function is defined as 

*(k,i/,/) = \F[s k At)}\ 2 

= ^o(KJ 2 + rU/27rit, 
1^t k m - ft)? + 1 



(5) 



i 



(6) 



the fully- anharmonic phonon frequency f^ v and lifetime 
Tk,i/ associated with a particular phonon mode (k, v) can 
be extracted by fitting to the Lorentzian peak shape of 
SED function \I/(k, ^, /). Here / = uj/2tt is the phonon 
frequency. 

On the other hand, the time derivative of the normal 
mode amplitude of mode (k, v) is given as 41 : 

3,n b ,N c , 

SuAt)= E J^«Sl'(t)e^(k,i/)ex P (-ik.r I ) (7) 

a,b,l r 

where a represents x,y,z directions, b is the index of ba- 
sis atoms, rib is the number of basis atoms in the chosen 
cell, / is the index of cells, N c is the total number of cells 
in the MD domain, and is the atomic mass of basis 
atom b. ify 1 is the a component of the velocity of basis 
atom b in cell I and r l is the equilibrium position of cell 
I. e^*(k, v) is the complex conjugate of the a component 
of the eigen-displacement of basis atom b associated with 
the ^th branch at a chosen wavevector k. Based on the 
time history of atomic velocities generated by MD sim- 
ulations and eigen-displacements from lattice dynamics 
(LD) calculations, the normal mode amplitude can be 
obtained according to Eq. ([7]). Then SED function can 
be constructed and fitted with Eq. (|6j) to extract f£ v 
and Tk,/v, as shown in Fig. [2j 

As mentioned in previous section, to minimize statisti- 
cal fluctuations due to MD simulations, five independent 
simulations are run for each case study. The arithmetic 
average of the corresponding SED functions are found to 
give better Lorentzian peak shapes. Therefore, all data in 
this work are extracted from the averaged SED functions. 
However, it should be pointed out that artificial broaden- 
ing of the peaks after averaging can be non-negligible if 
the peak positions from individual runs differ much. Such 
broadening will lead to underestimated phonon lifetime. 
To examine such effect, we assume two identical peaks are 
separated by a small amount 26. When they are super- 
positioned, the new peak centers at /. Therefore, the 
positions where half maximum of the new peak occurs 
are solutions to 
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(/-/-^) 2 + r 2 (/- / + £) 2 + r 2 (5 2 + r 2 
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FIG. 2: (Color online) Semilogarithmic plots of the sum of 
SED functions tf(k,/) = X^( k ^/) alon g r " K direction 
with k = k max /2. a) suspended SLG. b) supported SLG. The 
red dashed lines indicate phonon peak positions of suspended 
SLG. The numbers adjacent to each peak are lifetime values 
obtained through fitting to Eq. J6]). 



The new FWHM is then found as: 



2r' = 2rW2( f )2. 



5(^) 4 + 2(^) 2 + l. 



(9) 



In present study, it is usually observed that S/T < 0.4 
even for low frequency acoustic phonons. Therefore, the 
underestimation of phonon lifetime is expected to be well 
below 20%. This is confirmed by the agreement between 
the average of phonon lifetimes extracted from individual 
runs and that from averaged SED function. 

The thermal conductivity K\ can be obtained as a sum- 
mation of contributions from all phonon modes in the full 
Brillouin zone (BZ): 



^^^^^(k^Mk,!/), (io) 



where c v is the quantum per-mode volumetric specific 
heat of a harmonic oscillator, which decreases rapidly as 
frequency increases: 



ksx 2 exp(x) 
y[exp(x) - l] 2 



(11) 



Here V is the effective volume of the simulation domain, 
and x = Uuj/kBT. v g = duj/dk is the phonon group veloc- 
ity. The summation goes over all resolvable wavevector 
k's and all phonon branches v. 

One limitation of SED approach is that only a number 
of k points which comply with the periodicity of the MD 
domain can be resolved, leading to the finite resolution 
of the k-grid and cutoff of contribution from very long 
wavelength phonon modes ^ 39 i 42 . The exclusion of these 
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modes may lead to domain size effect of the thermal con- 
ductivity. Such effect generally presents in MD simula- 
tions with PBC specified rather than only exists in SED 
approach. However, due to the large amount of resolv- 
able phonon modes (k, v) in BZ associated with larger 
domains, study of domain size effect is not as feasible for 
SED approach. Nonetheless, a reasonably sized simula- 
tion domain should preserve the validness of at least the 
qualitative description of physical processes. 

One source of inaccuracy in our calculation may come 
from the fact that normal and Umklapp processes can 
both contribute to the broadened phonon peak in a non- 
distinguishable manner in MD simulations. While Umk- 
lapp process is directly resistive for thermal transport, 
normal process is only indirect - it increases the popu- 
lation of high wavevector phonons which can later en- 
counter Umklapp process and lead to thermal resistance. 
Compared to the simple addition of these two processes 
in our SED method, a more accurate treatment would 
be a consideration of these two processes separately and 
a solution of the full Boltzmann transport equatio n 15 ! 43 . 
However, that approach cannot take into account higher 
order phonon scattering processes or extrinsic scattering 
mechanisms such as substrate effects which are impor- 
tant for our problem. Considering the fact that a high 
wavevector phonon generated from a normal process will 
likely participate in the Umklapp process at a later time 
and contribute to thermal resistance, the direct inclusion 
of normal process in lifetime calculation should be ap- 
proximately valid. Also, the lifetime calculated from MD 
has already included the effects of non-equilbrium fluc- 
tuation of phonon distribution function, similar in spirit 
to that in the GK method 44 . 



III. RESULTS 
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A. Phonon lifetime 

The per-branch phonon lifetimes for both suspended 
and supported SLG at 304 K as a function of wavevec- 
tor k and frequency are shown in Fig. [3] (a) and (b), 
respectively. Because of symmetry, the phonon proper- 
ties in the full BZ can be reproduced by only studying 
the irreducible k-space, which is l/12 th of the 1 st BZ. 
However, because of the simplicity in discretization and 
specifying group velocity along x direction, we study the 
1 st quadrant of BZ instead. From Fig. [3] (a), we indeed 
find the distribution of lifetimes in suspended SLG to be 
reasonably symmetrical about high symmetry lines, con- 
sidering finite discretization and uncertainty in predicted 
lifetime values. This confirms the self-consistency of our 
approach. In addition, as seen in Fig [3] (a), the phonon 
dispersions in SLG with or without substrate do not differ 
much. Therefore, good correspondence between phonon 
modes in both cases can be well found to identify lifetime 
changes. 



FIG. 3: (Color online) a) Distribution of phonon lifetime in 
the first quadrant of 1 st BZ with the color representing life- 
times values. Also plotted are the phonon dispersions along 
T-K. b) Lifetime as a function of frequency for all phonon 
branches. The subscripts "s" and "p" are short for "sus- 
pended" and "supported", respectively. 



1. Phonons in suspended graphene 

Existing studies either suggest negligible— or 
dominant 1 - contribution from flexural phonons. A clear 
and correct picture would be crucial for understanding 
thermal transport in graphene. From Fig. [3j it is seen 
that, in suspended SLG, lifetime values of ZA phonons 
are in the range of sub 10-40 ps throughout BZ. The 
long lifetime of ZA phonons indicates relatively weak 
correlation between the in-plane and out-of-plane modes, 
which results in less scattering in ZA phonons^. It is 
interesting to see that out-of-plane optical (ZO) and 
ZA phonons have similarly large lifetimes, which can 
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FIG. 4: In-plane phonon lifetime in xy-SLG as compared to 
that in suspended SLG along Y — K direction. 



again be ascribed to the coupling strength between 
the in-plane and out-of-plane modes. The relatively 
weak dependence of ZA/ZO lifetime on the wavevector 
k is also an indication that their coupling to in-plane 
phonons are not very strong. As frequency increases, 
lifetime of ZA phonons first decreases then increases as 
k approaches BZ boundary. The relatively short lifetime 
of mid-frequency ZA phonons can be attributed to the 
fact that higher-order phonon scattering involving these 
modes is important and yields comparable scattering 
rates as that of three-phonon process, limiting the 
lifetime^ 5 -. 

To further explore the role of flexural phonons in sus- 
pended SLG, we carry out simulations on suspended 
SLG with vibrations restricted to within the x-y plane 
(xy-SLG) so that the atomic motion is completely 2D. 
In doing this, out-of-plane phonons (ZA and ZO) are 
completely removed from the system. Predictions from 
GK method (not shown) consistently suggest several 
times larger ni in xy-SLG. If only three-phonon scat- 
tering was considered, as suggested by Lindsay et al, 
ZA phonons would contribute about 80% of the thermal 
conductivity^ 5 -. Then the lifetime of remaining phonon 
modes should increase by at least four times to com- 
pensate the removal of phonon modes to result in an 
increased acj, according to Eq. ([TO]) . However, from 
SED analysis, as shown in the inset of Fig. [4j it is 
found that lifetime in xy-SLG is generally only two times 
longer than that in SLG in average. This implies the 
out-of-plane modes are coupled to in-plane modes with 
medium strength to leverage the scattering in in-plane 
and out-of-plane phonons in SLG. Such observation can 
only be understood when higher-order anharmonicities 
are taken into consideration. Due to the special symme- 
try of SLG, very stringent restrictions are imposed on 
the 3 rd order scattering processes involving out-of-plane 
phonon modes 15 . There are so few 3 rd order processes 
allowed that the higher-order scattering processes may 
yield considerable scattering rates and destroy the " pris- 
tine" transport of out-of-plane phonons^ 5 -, as supported 
by our simulation results. 

Away from the zone center, the overall lifetime is 
mainly around 20 ps and below in most regions of BZ for 
longitudinal acoustic (LA) and transverse acoustic (TA) 



phonons. As frequency increases, lifetime shortens due 
to more efficient Umklapp scattering. For LA and TA 
phonons near zone center with small frequencies, similar 
to ZA phonons, their lifetime tends to diverge as seen in 
Fig. [3](b). This is because these phonon modes have very 
small wavevectors that, due to the momentum and en- 
ergy conservation requirements, barely scatter with other 
phonons. 

As a side note, it is seen that lifetime in SLG is ac- 
tually not much higher than that in some solids with 
much lower thermal conductivity, such as solid argon 3 -^, 
silicon 33 and PbTe 36 . Such observation indicates that 
despite the unique 2D transport in SLG, phonons do not 
necessarily experience much less scattering. The large K\ 
of SLG is mainly due to the small atomic weight of car- 
bon atoms and the strong C-C bonds which lead to high 
phonon group velocity. 



2. Phonons in supported graphene 

When put on substrate, it is found that lifetimes are 
shortened throughout BZ and the full frequency range 
for all phonon branches. Namely, lifetimes of all acous- 
tic phonons are generally shortened to be less than 10 
ps. From the color map in Fig. [3] (a), it can be clearly 
seen that such reduction is strongest for phonon modes 
closer to zone center. Also the distribution of lifetime in 
BZ becomes less symmetrical. This indicates the intro- 
duction of substrate strongly destructs the long-range or- 
der of phonon transport through random scattering sites 
at the interface. From Fig. [3] (b) it is seen that the 
lifetime of low frequency phonons are mostly shortened. 
For frequency below 18 THz, the lifetime values mono- 
tonically decrease as frequency continuously reduces to 
zero, indicating the substrate couples most efficiently to 
low frequency acoustic phonons. One of the most signif- 
icant outcomes of such observation is that the thermal 
conductivity of supported graphene will not be sensitive 
to graphene flake size. This is because long wavelength 
near-zone-center phonons in larger samples will have min- 
imal lifetime due to substrate scattering and thus barely 
contribute to thermal transport. 

Recent Brillouin scattering experiment on multilayer 
graphene supported on SiC^/Si substrate 46 suggests the 
lifetime of near-zone-center LA/TA phonon is about 10- 
30 ps, which is shorter than that of suspended graphene, 
in line with our findings. In fact, due to the amor- 
phous nature of silica, the coordinations of surface sili- 
con/oxygen atoms cannot comply with the periodicity of 
SLG. As a result, the presence of the "irregularly" coor- 
dinated silicon/oxygen atoms at the interface will couple 
to the in-plane vibrations of carbon atoms as scattering 
centers and essentially hinder the in-plane phonons in 
SLG from transporting in a perfect 2D lattice. There- 
fore, to sustain long lifetime and thus high thermal con- 
ductivity of SLG, a substrate with better lattice match 
to SLG should be pursued. One model example is the 
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multiple-layer graphene. When all but one layers are 
treated as substrate, the effective "SLG" is supported on 
a " substrate" with perfectly- matched lattice that the LA 
and TA phonons are less affected^. Therefore a sup- 
ported double-layer graphene (DLG) or SLG supported 
on Boron Nitride (BN) substrate may turn out to be bet- 
ter for planar heat dissipation applications 47 . 

The lifetime reduction is especially strong for ZA 
modes. This is due to the fact that the presence of the 
substrate breaks various symmetries (mirror-reflection, 
translational, etc.) in SLG and alters the out-of-plane 
vibrations of SLG by introducing SLG-substrate scatter- 
ing which largely shortens the lifetime of ZA phonons. 
In this case, higher-order anharmonic interactions are no 
longer important. Similar arguments can also be made 
to explain the reduction of lifetime of ZO phonons. 



3. Longitudinal and transverse optical phonons 

The majority of longitudinal optical (LO) and trans- 
verse optical (TO) phonons have lifetimes about 1-4 ps at 
room temperature and reduced by about half when put 
on substrate. This is consistent with the optical phonon 
lifetimes of supported few-layer graphene measured by 
Raman spectroscopy 4 -"—. 



B. Phonon mean free path 
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As phonon spectral lifetimes are obtained for individ- 
ual phonon modes in the preceding section, it is straight- 
forward to compute phonon MFP based on its definition: 

7 duo / \ 

Iph = VgT = ~^T. (12) 

As seen in Fig. [3] (a), the phonon dispersion changes 
when SLG is put on substrate. This is due to the fact 
that the substrate induces internal stress and ripples in 
SLG, which leads to shift in phonon frequencies. The 
complete phonon dispersions before and after putting 
SLG on substrate are used to calculate the phonon group 
velocities, MFP's, thermal conductivities for suspended 
and supported SLG, respectively. MFP and per- mode 
thermal conductivity contribution of both suspended and 
supported SLG as a function of frequency are shown in 
Fig. [5l It is seen that MFP ranges from 10 nm to 600 
nm for all acoustic and ZO phonons while below 30 nm 
for LO and TO phonons. MFP's of all acoustic and ZO 
phonons are generally longer than that of LO and TO 
phonons' due to their longer lifetime and higher group 
velocity. MFP's of LA and TA phonons exhibit the ex- 
pected divergence when uo — » due to saturated group 
velocity and diverging lifetime. In contrast, no diver- 
gence is found for ZA modes because their group velocity 
v g ~ y/uJ appears to decrease faster than the divergence 
of lifetime as uo — >• 0. When put on substrate, the MFP 
values drop to less than 200 nm for all acoustic and ZO 



FIG. 5: (Color online) The frequency dependence of (a) 
MFP and (b) per-mode thermal conductivity for all phonon 
branches. 



phonons while those of LO and TO phonons are less af- 
fected. A comparison between the phonon dispersion re- 
lations before and after the presence of substrate doesn't 
suggest significant change in phonon group velocities, as 
seen in Fig. [3] (a). Therefore, according to Eq. (fT0]h 
the effect of substrate in limiting SLG thermal conduc- 
tivity should be mostly attributed to the introduction of 
SLG-substrate scattering that shortens lifetime. Due to 
the small MFP in supported SLG, again, we do not ex- 
pect graphene flake size to significantly affect the thermal 
conductivity as long as it is wider than 200 nm. 



C. Thermal conductivity decomposition 

From Fig. 0(b), it is seen that the majority of thermal 
conductivity contributions are from phonons with lower 
frequencies. This can be understood from Eq. ([T0|) that 
v g and r are large for lower frequencies and c v decreases 
with increasing frequency due to the change in phonon 
population. Therefore, despite their small group veloc- 
ity, Z A phonons still contribute about 1/4 of the thermal 
conductivity, which is not much less than that from LA 
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FIG. 6: (Color online) a) Contributions to thermal conduc- 
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contributions to thermal conductivity from ZA, TA, LA and 
ZO phonons. 



and TA phonons. Such contribution is mainly from mid- 
frequency phonons due to the different frequency depen- 
dence of ZA group velocity and lifetime. When put on 
substrate, the thermal conductivity contributed from all 
acoustic and ZO phonon branches are reduced, owing to 
the shortened lifetime due to SLG-substrate scattering 
mentioned in the above text. 

The relative importance of contribution to thermal 
conductivity from all phonon branches (except TO and 
LO) at elevated temperatures are computed using Eq. 
([TO]) and presented in Fig. [6] It is seen that all acoustic 
phonons contribute about 20-35% to the total thermal 
conductivity in suspended SLG. The relative contribu- 
tion from Z A phonons increases with decreasing tempera- 
ture as other phonons freeze out. When put on substrate, 
contribution from all branches to the total thermal con- 
ductivity reduces. It can also be seen that the reduction 
in ZA phonon contributions are most severe as its rel- 
ative contribution drops to be around 15%, comparable 
to that of ZO phonons. In contrast, the relative contri- 
butions from both LA and TA phonons are still above 
30% when supported. Therefore, the in-plane acoustic 
phonons are mainly responsible for the thermal transport 
in supported SLG, consistent with the theory of Seol et 
aj£. At higher temperature, K\ values associated with all 
acoustic phonons decrease because of stronger Umklapp 
scattering. At lower temperature, k\ values increase ow- 
ing to the activation of higher frequency phonon modes, 
which is only captured after QC. 

The predicted total thermal conductivity for sus- 
pended graphene before QC is 1626 W/m K at 300 
K, which agrees with k\ = 1695 W/m K obtained 
from Green-Kubo method. The experimentally measured 
thermal conductivity values for suspended SLG spread 
widely, which can be partly attributed to the varied 
graphene sample quality. As mentioned before, because 



FIG. 7: (Color online) Predicted total thermal conductivity 
values in comparison with experiments of suspended SLGf^, 
membrane 51 and SLG supported on Si02 substrate 6 . 

of the large uncertainty and spread data, no clear tem- 
perature trend or graphene flake size dependence can be 
demonstrated experimentally so far. Nonetheless, as seen 
in Fig. [7J the k\ values predicted from MD simulations 
are consistent with experiments for both suspended SLG 
and SLG on amorphous Si02 substrate. Such consistency 
and the successful prediction of thermal conductivity re- 
duction in supported SLG without adjustable parameters 
suggest the validity of the methodology and the results 
presented in this work. 



IV. FURTHER EXPLORATIONS 

A. Effects of interfacial bonding 

The above results for supported SLG are obtained un- 
der the assumption that the interfacial interactions be- 
tween SLG and substrate is of vdW type 20 . On the other 
hand, it is found that covalent bonds can form between 
graphene and SiC substrate^ 2 -. To explore the effects 
of interfacial bondings on the thermal transport in sup- 
ported SLG, we study a) lifetimes when LJ interactions 
are assumed to be 10 times stronger than used in previ- 
ous sections; b) lifetimes of graphene on silicon substrate 
with either vdW interactions or covalent bonds at the 
interface. For the latter case, we use LJ potentials with 
parameters ec-Si = 8.909 meV, crc-Si = 3.326 A to 
model vdW interactions while using OPT with mixing 
rules to model the C-Si covalent bonds. SED analysis 
is performed and the results are compared in T — K 
direction, as shown in Fig. [8] (a) and (b). As seen, 
when the vdW interactions are 10 times stronger or all 
bonds at the interface are covalent, the phonon scatter- 
ing in SLG due to substrate is much stronger, leading 
to significantly shortened lifetime and thus strongly re- 
duced thermal conductivity. According to the studies on 
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FIG. 8: (Color online) Lifetime of all phonon branches along 
r — K in supported SLG. a) SiCb substrate, when 10 times 
stronger LJ interactions assumed, b) Silicon substrate, when 
covalent interfacial interactions are assumed. 



charge density distribution at CNT/silicon interface^, 
the interactions between SLG and a substrate with non- 
passivated surface appear to be a mixture of non-bonded 
vdW interactions and covalent bonds. Therefore, in gen- 
eral, the value of thermal conductivity of supported SLG 
may be a weighted average of the values when both vdW 
and covalent bonds at the interface are assumed. 

On the other hand, we also find similar reduction in 
thermal conductivity when silicon substrate is used in the 
simulation with vdW interaction strength similar to that 
assumed in the preceding sections where Si02 substrate 
is used 37 . Therefore, we hypothesize that the interac- 
tion strength and geometry at SLG-substrate interface 
should be deterministic factors for thermal conductivity 
reduction in supported SLG. Our theory is supported by 
recent experiment that polymeric residue can bring the 
thermal conductivity of suspended DLG down to around 
600 W/m K 53 , which is similar to that found in SLG 
supported on Si02 substrate^. 



B. Effects of interatomic potentials 

It is known that REBO potentials 54 for hydrocar- 
bons tend to underestimate the thermal conductivity of 
graphene nanoribbons (GNR) and CNT's, possibly due 
to the inadequateness of these potentials in representing 
the anharmonicity in carbon-based systems. Nonethe- 
less, due to the popularity of REBO potentials in litera- 
ture, we also use them to perform MD simulations as well 
as SED analysis for both suspended and supported SLG 
then compare to the results from OPT, as summarized 
in Table HI As seen, despite the underestimation in total 
thermal conductivity values, REBO potentials can simi- 
larly reproduce relative importance of phonon branches 



TABLE I: Summary of thermal conductivity and its decom- 
position in suspended and supported SLG predicted by OPT 
and REBO potentials at 435 K. The units for ki and individ- 
ual mode contributions are (W/m K) and (%), respectively. 



Suspended 


ki ZA TA LA ZO 


OPT 


1082 23.6 31.6 30.0 13.9 


REBO 


509 33.8 27.7 30.5 7.9 


Supported 


ki ZA TA LA ZO 


OPT 


451 16.8 32.5 33.8 16.0 


REBO 


224 27.5 28.7 29.5 13.7 



in contributing to ki and the reduction of K\ when SLG 
is put on substrate. However, it does fail to predict sig- 
nificant reduction in ZA contribution in supported SLG. 
Therefore, the use of REBO potentials to study carbon- 
based systems should generally be valid in a qualitative 
sense while caution is advised. 



V. CONCLUSIONS 

In conclusion, we predict the phonon lifetime, MFP 
and thermal conductivity of both suspended and sup- 
ported SLG using MD simulations and SED analysis. 
Our results are found to be consistent with both opti- 
cal and thermal measurements for phonon lifetime and 
thermal conductivity of SLG in literature. In suspended 
SLG, out-of-plane phonons (ZA and ZO) are found to 
couple to the in-plane phonons in the Umklapp scatter- 
ing due to both 3 rd and higher-order anharmonic interac- 
tions, leading to around 25-30% contribution to thermal 
conductivity from ZA phonons. In the presence of sub- 
strate, the contributions to thermal conductivity from all 
acoustic and ZO phonon branches are reduced, owing to 
the SLG-substrate scattering and the breakdown of sym- 
metries in both in-plane and out-of-plane phonons. ZA 
phonon contributions to thermal conductivity are largely 
suppressed to be around 15% in supported SLG, making 
LA and TA phonons most responsible for thermal trans- 
port when SLG is put on amorphous Si02 substrate. We 
suggest to use substrates with closer lattice match to SLG 
for better planar heat dissipations. 
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Appendix A: Quantum Corrections 

In MD, quantum effects in specific heat and phonon 
occupation number are not captured due to its classi- 
cal nature. Therefore, MD simulations are generally not 
readily suitable to study thermal transport at tempera- 
tures much lower than the Debye temperature. While it 
is difficult to explicitly include quantum effects in MD 
simulations^- 5 , there are some quantum correction (QC) 
approaches available based on ad hoc physical arguments. 
According to them, the temperature correction is made 
by equating the total energies of the classical and quan- 
tum systems: 



Ec(T c )=E Q (T Q ), 



(Al) 



and the thermal conductivity correction is obtained as: 

dT c 



KQ = K C 



dT Q - 



(A2) 



While overall k values from such QCs may seem okey, the 
correctness of individual phonon mode properties can not 
be guaranteed 55 . Fortunately, we are able to obtain indi- 
vidual phonon properties from the spectral energy den- 
sity (SED) analysis, enabling QCs in a mode-dependent 
manner. The phonon lifetime is largely determined by 
the phonon distribution in the system. Therefore, to 
maximally ensure the validness of MD simulation in re- 
producing the mode-dependent phonon scattering, we 



perform QC by finding a quantum system at tempera- 
ture Tq that has the closest phonon distribution to the 
classical system at temperature Tc- That is to find a Tq 
that leads to the minimum of function 



M(T Q ,T C ) 



\fq(T Q ,w) - f c (T c ,u)\D(uj)duj 



(A3) 

for a choice of Tc- Here /q(Tq,o;) = l/[exp(xQ) — 1] and 
fc(Tc,uo) = 1/xc are the quantum and classical phonon 
occupation numbers, respectively, where Xi = huj/ksTi. 
D(u) is the phonon density of states that can be directly 
extracted from MD simulations through velocity- velocity 
autocorrelations. As an example, for suspended single- 
layer graphene (SLG) at Tc = 190 K, the equivalent 
quantum system is found to be at Tq = 304 K. 

Once Tq is found, we have tq(uo^ Vi Tq) w 
Tc(uo\t,v,Tc). By taking into account the quantum spe- 
cific heat 



k B XQ exp(xQ) 
F[exp(x Q ) - l] 2 



(A4) 



a correction to thermal conductivity can be obtained as: 



EE 



(A5) 
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